Analysis
Amanda Holzer, Arzie Bajrami, Alexandra Bagaïni, Rui Mata
2024-06-25
1 Scoping Review
This document visualizes the data for the scoping review.
1.1 Study Designs Across Domains
This plot shows the number of studies using a longitudinal or serial cross-sectional design across several domains.s
### plot domains with study design
### Function to create a single horizontal bar plot for counts of 1
create_combined_horizontal_bar <- function(df) {
columns_of_interest <- c('health', 'finance', 'political', 'crime', 'nature', 'nuclear', 'social')
# Create a long-format data frame for ggplot
df_long <- df %>%
tidyr::gather(key = "variable", value = "value", columns_of_interest) %>%
filter(value == 1)
# Arrange levels of 'variable' by count in ascending order
df_long$variable <- factor(df_long$variable, levels = names(sort(table(df_long$variable), decreasing = FALSE)))
# Define custom colors using hex color codes
custom_colors <- c("#440154", "#22A884") # Hex color codes without alpha
# Plotting the horizontal bar plot using ggplot2
ggplot(df_long, aes(x = value, y = variable, fill = study_design)) +
geom_bar(stat = "identity", position = position_stack(reverse = TRUE), color = NA) + # Reverse the stacking order
scale_fill_manual(values = setNames(custom_colors, unique(df_long$study_design))) + # Apply custom colors to study_design levels
theme_minimal() +
theme(panel.grid = element_blank(), panel.border = element_blank()) + # Remove grid lines and borders
labs(x = "Number of Studies", y = NULL, title = NULL) + # Remove y-axis label and title
xlab("Number of Studies") # Set x-axis label
}
# Apply the function to df_final
create_combined_horizontal_bar(df_final)
1.2 Measuring Risk Perception Across Domains
This plot shows the number of studies and how many times they measured risk perception across several domains.
###plot domain with measurement points-----
#create measurement category with 2, 3, and 3+
df_final <- df_final %>%
mutate(measure = ifelse(times_measured_1 == 2, "2",
ifelse(times_measured_1 == 3, "3", "3+")))
# Function to create a single horizontal bar plot for counts of 1
create_combined_horizontal_bar <- function(df) {
columns_of_interest <- c('health', 'finance', 'political', 'crime', 'nature', 'nuclear', 'social')
# Create a long-format data frame for ggplot
df_long <- df %>%
gather(key = "variable", value = "value", columns_of_interest) %>%
filter(value == 1) # Filter only rows where the value is 1
# Arrange levels of 'variable' by count in ascending order
df_long$variable <- factor(df_long$variable, levels = names(sort(table(df_long$variable), decreasing = FALSE)))
# Plotting the horizontal bar plot using ggplot2 with Viridis color scale
ggplot(df_long, aes(x = value, y = variable, fill = measure)) +
geom_bar(stat = "identity", position = position_stack(reverse = TRUE), color = NA) + # Reverse the stacking order
scale_fill_viridis_d() + # Use Viridis color palette for measure variable
theme_minimal() +
theme(panel.grid = element_blank(), panel.border = element_blank()) + # Remove grid lines and borders
labs(x = "Number of Studies", y = NULL, title = NULL) + # Remove y-axis label and title
xlab("Number of Studies") # Set x-axis label
}
# Apply the function to df_final
create_combined_horizontal_bar(df_final)
1.3 Risk Formulation Across Domains
This plot illustrates the different risk formulations used to measure risk perception across various domains.
#kind of risk (risk, concern, worry)
create_combined_horizontal_bar <- function(df) {
columns_of_interest <- c('health', 'finance', 'political', 'crime', 'nature', 'nuclear', 'social')
# Create a long-format data frame for ggplot
df_long <- df %>%
gather(key = "variable", value = "value", columns_of_interest) %>%
filter(value == 1) # Filter only rows where the value is 1
# Plotting the horizontal bar plot using ggplot2 with Viridis color scale
ggplot(df_long, aes(x = variable, fill = as.factor(risk_1))) +
geom_bar(position = "fill", stat = "count", color = "white") +
coord_flip() + # Flip the coordinates for a horizontal bar plot
theme_minimal() +
theme(panel.grid = element_blank(), panel.border = element_blank(),
axis.text.x = element_blank(), axis.title.x = element_blank()) + # Remove x-axis labels and title
labs(x = NULL, y = NULL, title = NULL) + # Remove axis labels and change title
scale_fill_viridis_d(name = "Risk Formulation") # Use Viridis color palette and customize legend title
}
# Apply the function to df_final
create_combined_horizontal_bar(df_final)
1.4 Event Occurence Across Domains
This plot shows the occurrence of events across various domains.
# Define the function
create_combined_horizontal_bar <- function(df) {
columns_of_interest <- c('health', 'finance', 'political', 'crime', 'nature', 'nuclear', 'social')
# Create a long-format data frame for ggplot
df_long <- df %>%
gather(key = "variable", value = "value", columns_of_interest) %>%
filter(value == 1) # Filter only rows where the value is 1
# Create a new factor combining intervention and exposure columns into "event" and "no event"
df_long$event_category <- ifelse(
df_long$intervention_yesno_1 == 0 & df_long$exposure_yesno_1 == 0,
"No Event",
"Event"
)
# Plotting the horizontal bar plot using ggplot2 with Viridis color scale
ggplot(df_long, aes(x = variable, fill = event_category)) +
geom_bar(position = "fill", stat = "count", color = "white") +
coord_flip() + # Flip the coordinates for a horizontal bar plot
theme_minimal() +
theme(panel.grid = element_blank(), panel.border = element_blank(),
axis.text.x = element_blank(), axis.title.x = element_blank()) + # Remove x-axis labels and title
labs(x = NULL, y = NULL, title = NULL) + # Remove axis labels and change title
scale_fill_viridis_d(name = "Event Occurence", labels = c("No Event", "Event")) # Use Viridis color palette and customize legend title
}
# Apply the function to df_final
create_combined_horizontal_bar(df_final)
1.5 Demographics Table Placeholder
Participant characteristics and demographics.
#Table2 Temporal Stability
2.1 Model Fitting and Evaluation - Meta-Analytical Stability and Change Model (MASC)
2.1.1 Model 1: Fitting with Correlation, Interval, and Size
# Define family
family <- brmsfamily(
family = "student",
link = "identity"
)
# Define the formula
formula_m1 <- bf(
cor_val | resp_se(sei, sigma = TRUE) ~ rel * (change * ((stabch^interval_val) - 1) + 1),
nlf(rel ~ inv_logit(logitrel)),
nlf(change ~ inv_logit(logitchange)),
nlf(stabch ~ inv_logit(logitstabch)),
logitrel ~ 1,
logitchange ~ 1,
logitstabch ~ 1,
nl = TRUE
)
# Define the weakly informative priors
priors <-
prior(normal(0,1), nlpar="logitrel", class = "b") +
prior(normal(0,1), nlpar="logitchange", class = "b") +
prior(normal(0,1), nlpar="logitstabch", class = "b") +
prior(cauchy(0,1), class = "sigma")
# Fit the model
fit_masc_m1 <- brm(
formula = formula_m1,
prior = priors,
family = family,
data = df,
cores = 2,
chains = 2,
iter = 6000,
warmup = 2000,
control = list(max_treedepth = 10, adapt_delta = 0.95),
seed = 1299
)## Compiling Stan program...## Start sampling# Model summary
summary(fit_masc_m1)## Family: student
## Links: mu = identity; sigma = identity; nu = identity
## Formula: cor_val | resp_se(sei, sigma = TRUE) ~ rel * (change * ((stabch^interval_val) - 1) + 1)
## rel ~ inv_logit(logitrel)
## change ~ inv_logit(logitchange)
## stabch ~ inv_logit(logitstabch)
## logitrel ~ 1
## logitchange ~ 1
## logitstabch ~ 1
## Data: df (Number of observations: 101)
## Draws: 2 chains, each with iter = 6000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## logitrel_Intercept 0.67 0.17 0.38 1.04 1.00 3867
## logitchange_Intercept -0.76 0.58 -1.84 0.58 1.00 3667
## logitstabch_Intercept -0.36 1.09 -2.44 1.75 1.00 3479
## Tail_ESS
## logitrel_Intercept 4310
## logitchange_Intercept 2581
## logitstabch_Intercept 3832
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.19 0.02 0.15 0.23 1.00 4261 4046
## nu 17.66 12.01 4.20 50.21 1.00 4216 4488
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Plot conditional effects
plot(conditional_effects(fit_masc_m1), points = TRUE)
# Trace plots and parameter estimates
plot(fit_masc_m1, N = 5, ask = TRUE)
# Posterior predictive checks
pp_check(fit_masc_m1, type = "dens_overlay", ndraws = 100)
pp_check(fit_masc_m1, type = "stat", stat = "mean", ndraws = 1000, binwidth = .001)
# LOO and diagnostics
model_loo_m1 <- loo(fit_masc_m1, save_psis = TRUE, cores = 2)
plot(model_loo_m1, diagnostic = "k")
plot(model_loo_m1, diagnostic = "n_eff")
2.1.2 Model 2a: Fitting with Demographics - Age, Age^2, Gender
# Define family
family <- brmsfamily(
family = "student",
link = "identity"
)
# Define the formula
formula_m2a <- bf(
cor_val | resp_se(sei, sigma = TRUE) ~ rel * (change * ((stabch^interval_val) - 1) + 1),
nlf(rel ~ inv_logit(logitrel)),
nlf(change ~ inv_logit(logitchange)),
nlf(stabch ~ inv_logit(logitstabch)),
logitrel ~ 1 + age_dec_c + age_dec_c2 + female_c,
logitchange ~ 1 + age_dec_c + age_dec_c2 + female_c,
logitstabch ~ 1 + age_dec_c + age_dec_c2 + female_c,
nl = TRUE
)
# Define the weakly informative priors
priors <-
prior(normal(0,1), nlpar="logitrel", class = "b") +
prior(normal(0,1), nlpar="logitchange", class = "b") +
prior(normal(0,1), nlpar="logitstabch", class = "b") +
prior(cauchy(0,1), class = "sigma")
# Fit the model
fit_masc_m2a <- brm(
formula = formula_m2a,
prior = priors,
family = family,
data = df,
cores = 2,
chains = 2,
iter = 6000,
warmup = 2000,
control = list(max_treedepth = 10, adapt_delta = 0.95),
seed = 1299
)## Compiling Stan program...## Start sampling# Model summary
summary(fit_masc_m2a)## Family: student
## Links: mu = identity; sigma = identity; nu = identity
## Formula: cor_val | resp_se(sei, sigma = TRUE) ~ rel * (change * ((stabch^interval_val) - 1) + 1)
## rel ~ inv_logit(logitrel)
## change ~ inv_logit(logitchange)
## stabch ~ inv_logit(logitstabch)
## logitrel ~ 1 + age_dec_c + age_dec_c2 + female_c
## logitchange ~ 1 + age_dec_c + age_dec_c2 + female_c
## logitstabch ~ 1 + age_dec_c + age_dec_c2 + female_c
## Data: df (Number of observations: 101)
## Draws: 2 chains, each with iter = 6000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## logitrel_Intercept 0.49 0.17 0.16 0.85 1.00 5441
## logitrel_age_dec_c -0.00 0.26 -0.47 0.60 1.00 2283
## logitrel_age_dec_c2 0.22 0.24 -0.00 0.90 1.00 1632
## logitrel_female_c 0.49 0.44 -0.36 1.35 1.00 7709
## logitchange_Intercept -0.74 0.80 -2.29 0.86 1.00 3341
## logitchange_age_dec_c -0.33 0.59 -1.65 0.92 1.00 2786
## logitchange_age_dec_c2 -0.03 0.33 -0.61 0.71 1.00 2401
## logitchange_female_c 0.16 0.80 -1.43 1.70 1.00 6218
## logitstabch_Intercept 0.60 1.18 -1.94 2.62 1.00 4362
## logitstabch_age_dec_c 0.16 0.86 -1.57 1.82 1.00 5099
## logitstabch_age_dec_c2 -0.77 0.65 -2.00 0.59 1.00 2338
## logitstabch_female_c 0.02 0.99 -1.94 1.91 1.00 8253
## Tail_ESS
## logitrel_Intercept 4035
## logitrel_age_dec_c 2287
## logitrel_age_dec_c2 1000
## logitrel_female_c 5304
## logitchange_Intercept 4888
## logitchange_age_dec_c 2657
## logitchange_age_dec_c2 2167
## logitchange_female_c 4874
## logitstabch_Intercept 5590
## logitstabch_age_dec_c 5164
## logitstabch_age_dec_c2 2477
## logitstabch_female_c 5755
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.18 0.02 0.13 0.23 1.00 3366 3896
## nu 13.96 11.18 2.85 44.29 1.00 3860 4131
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Plot conditional effects
plot(conditional_effects(fit_masc_m2a), points = TRUE)
# Trace plots and parameter estimates
plot(fit_masc_m2a, N = 5, ask = TRUE)
# Posterior predictive checks
pp_check(fit_masc_m2a, type = "dens_overlay", ndraws = 100)
pp_check(fit_masc_m2a, type = "stat", stat = "mean", ndraws = 1000, binwidth = .001)
# Leave-One-Out (LOO) and diagnostics
model_loo_m2a <- loo(fit_masc_m2a, save_psis = TRUE, cores = 2)
plot(model_loo_m2a, diagnostic = "k")
plot(model_loo_m2a, diagnostic = "n_eff")
2.1.3 Model 2b: Fitting with Measurement Characteristics - Domain, Item, Event
# Define family
family <- brmsfamily(
family = "student",
link = "identity"
)
# Define the formula
formula_m2b <- bf(
cor_val | resp_se(sei, sigma = TRUE) ~ rel * (change * ((stabch^interval_val) - 1) + 1),
nlf(rel ~ inv_logit(logitrel)),
nlf(change ~ inv_logit(logitchange)),
nlf(stabch ~ inv_logit(logitstabch)),
logitrel ~ 1 + domain + item + event,
logitchange ~ 1 + domain + event,
logitstabch ~ 1 + domain + event,
nl = TRUE
)
# Define the weakly informative priors
priors <-
prior(normal(0,1), nlpar="logitrel", class = "b") +
prior(normal(0,1), nlpar="logitchange", class = "b") +
prior(normal(0,1), nlpar="logitstabch", class = "b") +
prior(cauchy(0,1), class = "sigma")
# Fit the model
fit_masc_m2b <- brm(
formula = formula_m2b,
prior = priors,
family = family,
data = df,
cores = 2,
chains = 2,
iter = 6000,
warmup = 2000,
control = list(max_treedepth = 10, adapt_delta = 0.95),
seed = 1299
)## Compiling Stan program...## Start sampling# Model summary
summary(fit_masc_m2b)## Family: student
## Links: mu = identity; sigma = identity; nu = identity
## Formula: cor_val | resp_se(sei, sigma = TRUE) ~ rel * (change * ((stabch^interval_val) - 1) + 1)
## rel ~ inv_logit(logitrel)
## change ~ inv_logit(logitchange)
## stabch ~ inv_logit(logitstabch)
## logitrel ~ 1 + domain + item + event
## logitchange ~ 1 + domain + event
## logitstabch ~ 1 + domain + event
## Data: df (Number of observations: 101)
## Draws: 2 chains, each with iter = 6000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## logitrel_Intercept 1.22 0.37 0.61 2.06 1.00 3725
## logitrel_domainother -0.14 0.14 -0.40 0.13 1.00 5139
## logitrel_itemsingle 0.43 0.32 -0.07 1.19 1.00 4028
## logitrel_eventevent -0.68 0.16 -1.01 -0.40 1.00 5017
## logitchange_Intercept -0.73 0.76 -2.12 0.83 1.00 3907
## logitchange_domainother -0.43 0.76 -1.83 1.22 1.00 3093
## logitchange_eventevent -0.31 0.65 -1.58 1.26 1.00 2943
## logitstabch_Intercept 0.04 0.93 -1.84 1.81 1.00 5493
## logitstabch_domainother 0.63 1.00 -1.39 2.39 1.00 3501
## logitstabch_eventevent 0.26 0.95 -1.61 2.10 1.00 3289
## Tail_ESS
## logitrel_Intercept 2885
## logitrel_domainother 4837
## logitrel_itemsingle 3159
## logitrel_eventevent 4704
## logitchange_Intercept 5090
## logitchange_domainother 3986
## logitchange_eventevent 2863
## logitstabch_Intercept 5217
## logitstabch_domainother 4837
## logitstabch_eventevent 5142
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.14 0.02 0.09 0.19 1.00 3524 4407
## nu 10.25 9.67 2.26 36.80 1.00 3587 4955
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Plot conditional effects
plot(conditional_effects(fit_masc_m2b), points = TRUE)
# Trace plots and parameter estimates
plot(fit_masc_m2b, N = 5, ask = TRUE)
# Posterior predictive checks
pp_check(fit_masc_m2b, type = "dens_overlay", ndraws = 100)
pp_check(fit_masc_m2b, type = "stat", stat = "mean", ndraws = 1000, binwidth = .001)
# Leave-One-Out (LOO) and diagnostics
model_loo_m2b <- loo(fit_masc_m2b, save_psis = TRUE, cores = 2)
plot(model_loo_m2b, diagnostic = "k")
plot(model_loo_m2b, diagnostic = "n_eff")
2.1.4 Model 3: Fitting with Demographics and Measurement Characteristics - Domain, Item, Event, Formulation
# Define the family
family <- brmsfamily(
family = "student",
link = "identity"
)
# Define the formula
formula_m3 <- bf(
cor_val | resp_se(sei, sigma = TRUE) ~ rel * (change * ((stabch^interval_val) - 1) + 1),
nlf(rel ~ inv_logit(logitrel)),
nlf(change ~ inv_logit(logitchange)),
nlf(stabch ~ inv_logit(logitstabch)),
logitrel ~ 1 + age_dec_c + age_dec_c2 + female_c + domain + item + event + formulation,
logitchange ~ 1 + age_dec_c + age_dec_c2 + female_c + domain + item + event + formulation,
logitstabch ~ 1 + age_dec_c + age_dec_c2 + female_c + domain + item + event + formulation,
nl = TRUE
)
# Define the weakly informative priors
priors <-
prior(normal(0, 1), nlpar = "logitrel", class = "b") +
prior(normal(0, 1), nlpar = "logitchange", class = "b") +
prior(normal(0, 1), nlpar = "logitstabch", class = "b") +
prior(cauchy(0, 1), class = "sigma")
# Fit the model
fit_masc_m3 <- brm(
formula = formula_m3,
prior = priors,
family = family,
data = df,
cores = 2,
chains = 2,
iter = 6000,
warmup = 2000,
control = list(max_treedepth = 10, adapt_delta = 0.95),
seed = 1299
)## Compiling Stan program...## Start sampling# MODEL 3 EVAL: MCMC DIAGNOSTICS --------------------------------------------------------
# Model summary
summary(fit_masc_m3)## Family: student
## Links: mu = identity; sigma = identity; nu = identity
## Formula: cor_val | resp_se(sei, sigma = TRUE) ~ rel * (change * ((stabch^interval_val) - 1) + 1)
## rel ~ inv_logit(logitrel)
## change ~ inv_logit(logitchange)
## stabch ~ inv_logit(logitstabch)
## logitrel ~ 1 + age_dec_c + age_dec_c2 + female_c + domain + item + event + formulation
## logitchange ~ 1 + age_dec_c + age_dec_c2 + female_c + domain + item + event + formulation
## logitstabch ~ 1 + age_dec_c + age_dec_c2 + female_c + domain + item + event + formulation
## Data: df (Number of observations: 101)
## Draws: 2 chains, each with iter = 6000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## logitrel_Intercept 1.27 0.36 0.67 2.09 1.00
## logitrel_age_dec_c 0.06 0.15 -0.18 0.41 1.00
## logitrel_age_dec_c2 0.01 0.08 -0.11 0.19 1.00
## logitrel_female_c 0.24 0.42 -0.57 1.06 1.00
## logitrel_domainother -0.01 0.15 -0.27 0.30 1.00
## logitrel_itemsingle 0.17 0.29 -0.33 0.81 1.00
## logitrel_eventevent -0.42 0.15 -0.74 -0.14 1.00
## logitrel_formulationconcern 0.13 0.63 -1.01 1.50 1.00
## logitrel_formulationthreat -0.32 0.35 -0.99 0.40 1.00
## logitrel_formulationworry 0.80 0.37 0.11 1.59 1.00
## logitchange_Intercept -0.28 0.89 -1.96 1.50 1.00
## logitchange_age_dec_c 0.38 0.68 -0.92 1.80 1.00
## logitchange_age_dec_c2 0.49 0.61 -0.68 1.79 1.00
## logitchange_female_c 0.06 0.94 -1.85 1.90 1.00
## logitchange_domainother -0.08 0.78 -1.66 1.47 1.00
## logitchange_itemsingle 0.06 0.95 -1.84 1.85 1.00
## logitchange_eventevent 0.24 0.76 -1.35 1.70 1.00
## logitchange_formulationconcern 0.25 0.94 -1.58 2.13 1.00
## logitchange_formulationthreat 0.57 0.84 -1.16 2.20 1.00
## logitchange_formulationworry 0.30 0.79 -1.24 1.90 1.00
## logitstabch_Intercept 0.04 0.89 -1.73 1.77 1.00
## logitstabch_age_dec_c 0.16 0.61 -1.31 1.19 1.00
## logitstabch_age_dec_c2 0.02 0.47 -1.14 0.73 1.00
## logitstabch_female_c -0.13 0.91 -1.92 1.73 1.00
## logitstabch_domainother 0.14 0.79 -1.42 1.68 1.00
## logitstabch_itemsingle 0.53 0.91 -1.41 2.18 1.00
## logitstabch_eventevent 0.25 0.78 -1.31 1.74 1.00
## logitstabch_formulationconcern -0.22 0.90 -1.98 1.53 1.00
## logitstabch_formulationthreat -0.59 0.89 -2.31 1.22 1.00
## logitstabch_formulationworry -0.52 0.77 -2.05 1.03 1.00
## Bulk_ESS Tail_ESS
## logitrel_Intercept 3095 3052
## logitrel_age_dec_c 2569 1699
## logitrel_age_dec_c2 1977 1351
## logitrel_female_c 6872 5173
## logitrel_domainother 3286 4718
## logitrel_itemsingle 2673 2344
## logitrel_eventevent 3268 4174
## logitrel_formulationconcern 3785 4171
## logitrel_formulationthreat 4152 4037
## logitrel_formulationworry 2849 4353
## logitchange_Intercept 4108 5699
## logitchange_age_dec_c 2066 4377
## logitchange_age_dec_c2 1266 2138
## logitchange_female_c 8640 5622
## logitchange_domainother 3816 4717
## logitchange_itemsingle 2177 4058
## logitchange_eventevent 2106 3880
## logitchange_formulationconcern 7220 5937
## logitchange_formulationthreat 6172 5344
## logitchange_formulationworry 4792 4735
## logitstabch_Intercept 3434 4990
## logitstabch_age_dec_c 1762 3091
## logitstabch_age_dec_c2 1165 1901
## logitstabch_female_c 8073 6119
## logitstabch_domainother 4316 5266
## logitstabch_itemsingle 1904 3466
## logitstabch_eventevent 2048 3512
## logitstabch_formulationconcern 5528 5990
## logitstabch_formulationthreat 5513 5518
## logitstabch_formulationworry 3136 4218
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.10 0.03 0.05 0.16 1.00 1675 3489
## nu 4.76 5.30 1.40 19.48 1.00 2101 3951
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Plot conditional effects
plot(conditional_effects(fit_masc_m3), points = T)
# Trace plots & parameter estimates
plot(fit_masc_m3, N = 5, ask = TRUE)
# MODEL 3 EVAL: PP CHECKS --------------------------------------------------------
pp_check(fit_masc_m3, type = "dens_overlay", ndraws = 100)
pp_check(fit_masc_m3, type = "stat", stat = "mean", ndraws = 1000, binwidth = .001)
pp_check(fit_masc_m3, type = "stat", stat = "sd", ndraws = 500, binwidth = .001)
pp_check(fit_masc_m3, type = "stat", stat = "median", ndraws = 500, binwidth = .001)
pp_check(fit_masc_m3, type = "stat", stat = "mad", ndraws = 500, binwidth = .001)
pp_check(fit_masc_m3, type = "stat_2d")## Using all posterior draws for ppc type 'stat_2d' by default.
pp_check(fit_masc_m3, type = "scatter_avg")## Using all posterior draws for ppc type 'scatter_avg' by default.
# MODEL 3 EVAL: LOO --------------------------------------------------------
# LOO & Pareto K
model_loo_m3 <- loo(fit_masc_m3, save_psis = TRUE, cores = 2)
plot(model_loo_m3, diagnostic = "k")
plot(model_loo_m3, diagnostic = "n_eff")
# LOO PIT
w <- weights(model_loo_m3$psis_object)
ppc_loo_pit_overlay(
y = fit_masc_m3$data$cor_val,
yrep = posterior_predict(fit_masc_m3),
lw = w
)## NOTE: The kernel density estimate assumes continuous observations and is not optimal for discrete observations.
ppc_loo_pit_qq(
y = fit_masc_m3$data$cor_val,
yrep = posterior_predict(fit_masc_m3),
lw = w
)
# MODEL EVAL: LOO - COMPARISON OF ALL MODELS ----------------------------------------------------
loo1 <- loo(fit_masc_m1)
loo2a <- loo(fit_masc_m2a)
loo2b <- loo(fit_masc_m2b)
loo3 <- loo(fit_masc_m3)
loo_compare(loo1, loo2a, loo2b, loo3)## elpd_diff se_diff
## fit_masc_m3 0.0 0.0
## fit_masc_m2b -7.1 5.3
## fit_masc_m1 -21.3 6.1
## fit_masc_m2a -21.7 5.92.2 Predictions Time
2.2.1 Plot correlation and interval
nd <- crossing(domain= NA,
sei = 0.1,
item= NA,
event=NA,
formulation=NA,
age_dec_c= 0,
age_dec_c2= 0,
female_c= 0,
interval_val=seq(0,5, by = .1))
epred_draws_df <- nd %>%
add_epred_draws(fit_masc_m3, re_formula = NA)
ggplot(epred_draws_df) +
stat_lineribbon(alpha = 1/4, point_interval = "mean_hdci", aes(x = interval_val, y = .epred)) +
geom_point(data= df, aes(x=interval_val, y= cor_val), size = 1.5, color = "grey20", alpha = 0.7) +
theme_minimal() +
labs(y = "Retest Correlation", x = "Retest Interval (Years)", color = "", linetype = "", fill = "",
title = "Behaviour") +
theme(strip.placement = "outside",
legend.position = "none",
legend.margin = margin(-.5,0,0,0, unit="cm"),
legend.spacing.y = unit(0.15, 'cm'),
legend.key.width = unit(1, "cm"),
legend.key.size = unit(.3, "cm"),
legend.text = element_text(size = 9, color = "grey20"),
text = element_text(size = 9, color = "grey40"),
axis.text.y = element_text(vjust=seq(0,1, length.out = 5)),
axis.text.x = element_text(hjust=c(0,1)),
title = element_text(size = 9, color = "grey20"),
plot.tag = element_text(size = 11, face = "bold", color = "grey20"),
panel.spacing = unit(.5, "lines"),
plot.title = element_blank(),
panel.grid = element_blank(),
plot.title.position = "plot",
plot.margin = margin(b = 5, r = 5, l = 5),
panel.background = element_rect(color = "grey75", fill = NA, size = .4)) +
scale_fill_manual(values = c("#502a7b","#502a7c","#502a7a")) +
guides(color = guide_legend(override.aes = list(size = .75)),
fill = guide_legend(override.aes = list(size = .75)),
size = "none", linetype = guide_legend(override.aes = list(size = .75))) +
coord_cartesian(ylim = c(0, 1), xlim = c(0,5)) +
scale_y_continuous(breaks = seq(0,1,0.25)) +
scale_x_continuous(breaks = c(0,1,2,3,4,5))
2.2.2 Plot correlation and interval by event
nd <- crossing(domain= NA,
sei = 0.1,
item= NA,
event= unique(df$event),
formulation= NA,
age_dec_c= 0,
age_dec_c2= 0,
female_c= 0,
interval_val=seq(0,5, by = .1))
epred_draws_df <- nd %>%
add_epred_draws(fit_masc_m3, re_formula = NA)
# epred_draws_dom <- epred_draws_df %>%
# group_by(interval_val, health_subdomain) %>%
# mean_hdci(.epred,.width = c(.95,.8,.5)) %>%
# pivot_wider(names_from = .width, values_from = c(.lower,.upper))
# Set order and capitalise first letter
epred_draws_df$event <- str_to_title(epred_draws_df$event)
epred_draws_df$event<- factor(epred_draws_df$event, levels = c("None","Event"))
df_plot <- df %>%
mutate(event = str_to_title(event))
df_plot$event<- factor(df_plot$event, levels = c("None","Event"))
# Plot
p_event <- ggplot(epred_draws_df) +
stat_lineribbon(alpha = 1/4, point_interval = "mean_hdci", aes(x = interval_val, y = .epred)) +
geom_point(data= df_plot, aes(x=interval_val, y= cor_val), size = 1.5, color = "grey20", alpha = 0.7) +
# geom_line(data = epred_draws_agg,
# aes(x = time_diff_dec*10, y = .epred),
# color = "grey95",
# size = .5) +
# geom_line(data = epred_draws_dom,
# aes(x = time_diff_dec*10, y = .epred, linetype = domain_name),
# color = "#e07f00",
# linewidth = .25) +
# geom_text_repel(data = lbl_dot_df,
# aes(x = time_diff_dec*10, y = .epred, label = label),
# family = "Source Sans 3", size = 2.5,
# min.segment.length = 0,
# segment.color = "grey50",
# segment.size = .25,
# box.padding = 0.5,
# nudge_x = .5,
# nudge_y = c(.1, -.05)
# ) +
facet_wrap(.~event) +
theme_minimal() +
labs(y = "Retest Correlation", x = "Retest Interval (Years)", color = "", linetype = "", fill = "",
title = "Behaviour") +
theme(strip.placement = "outside",
legend.position = "none", # c(0,0) bottom left, c(1,1) top-right.
legend.margin = margin(-.5,0,0,0, unit="cm"),
legend.spacing.y = unit(0.15, 'cm'),
legend.key.width = unit(1, "cm"),
legend.key.size = unit(.3, "cm"),
# plot.tag.position = c(0,.8),
legend.text = element_text(size = 10, color = "grey20"),
text = element_text(size = 12, color = "grey40"),
axis.text.y = element_text( vjust=seq(0,1, length.out = 5)),
axis.text.x = element_text( hjust=c(0,1)),
title = element_text(size = 12, color = "grey20"),
plot.tag = element_text( size = 12, face = "bold", color = "grey20"),
panel.spacing = unit(.5, "lines"),
plot.title = element_blank(),
panel.grid = element_blank(),
plot.title.position = "plot",
plot.margin = margin(b = 5, r = 5, l = 5),
panel.background = element_rect(color = "grey75", fill = NA, size = .4),
strip.text = element_text(face = "bold"),
strip.text.y = element_text(size = 12),
strip.text.x = element_text(size = 12)) +
scale_fill_manual(values = c("#502a7b","#502a7c","#502a7a" ))+
guides(color = guide_legend(override.aes = list(size = .75)),
fill = guide_legend(override.aes = list(size = .75)),
size = "none", linetype = guide_legend(override.aes = list(size = .75))) +
coord_cartesian(ylim = c(0, 1), xlim = c(0,5))+
scale_y_continuous(breaks = seq(0,1,0.25)) +
scale_x_continuous(breaks = c(0,1,2,3,4,5))
p_event
ggsave(plot = p_event,
filename = here::here("analysis", "masc_pred_event.jpeg"),
dpi = 300, width = 28, height = 30, units = "cm")2.2.3 Plot correlation and interval by domain
nd <- crossing(domain= unique(df$domain),
sei = 0.1,
item= NA,
event= NA,
formulation= NA,
age_dec_c= 0,
age_dec_c2= 0,
female_c= 0,
interval_val=seq(0,5, by = .1))
epred_draws_df <- nd %>%
add_epred_draws(fit_masc_m3, re_formula = NA)
# epred_draws_dom <- epred_draws_df %>%
# group_by(interval_val, health_subdomain) %>%
# mean_hdci(.epred,.width = c(.95,.8,.5)) %>%
# pivot_wider(names_from = .width, values_from = c(.lower,.upper))
ggplot(epred_draws_df) +
stat_lineribbon(alpha = 1/4, point_interval = "mean_hdci", aes(x = interval_val, y = .epred)) +
geom_point(data= df, aes(x=interval_val, y= cor_val), size = 1.5, color = "grey20", alpha = 0.7) +
# geom_line(data = epred_draws_agg,
# aes(x = time_diff_dec*10, y = .epred),
# color = "grey95",
# size = .5) +
# geom_line(data = epred_draws_dom,
# aes(x = time_diff_dec*10, y = .epred, linetype = domain_name),
# color = "#e07f00",
# linewidth = .25) +
# geom_text_repel(data = lbl_dot_df,
# aes(x = time_diff_dec*10, y = .epred, label = label),
# family = "Source Sans 3", size = 2.5,
# min.segment.length = 0,
# segment.color = "grey50",
# segment.size = .25,
# box.padding = 0.5,
# nudge_x = .5,
# nudge_y = c(.1, -.05)
# ) +
facet_wrap(.~str_to_title(domain)) +
theme_minimal() +
labs(y = "Retest Correlation", x = "Retest Interval (Years)", color = "", linetype = "", fill = "",
title = "Behaviour") +
theme(strip.placement = "outside",
legend.position = "none", # c(0,0) bottom left, c(1,1) top-right.
legend.margin = margin(-.5,0,0,0, unit="cm"),
legend.spacing.y = unit(0.15, 'cm'),
legend.key.width = unit(1, "cm"),
legend.key.size = unit(.3, "cm"),
# plot.tag.position = c(0,.8),
legend.text = element_text(size = 10, color = "grey20"),
text = element_text(size = 12, color = "grey40"),
axis.text.y = element_text( vjust=seq(0,1, length.out = 5)),
axis.text.x = element_text( hjust=c(0,1)),
title = element_text(size = 12, color = "grey20"),
plot.tag = element_text( size = 12, face = "bold", color = "grey20"),
panel.spacing = unit(.5, "lines"),
plot.title = element_blank(),
panel.grid = element_blank(),
plot.title.position = "plot",
plot.margin = margin(b = 5, r = 5, l = 5),
panel.background = element_rect(color = "grey75", fill = NA, size = .4),
strip.text = element_text(face = "bold"),
strip.text.y = element_text(size = 12),
strip.text.x = element_text(size = 12)) +
scale_fill_manual(values = c("#502a7b","#502a7c" , "#502a7a" )) +
guides(color = guide_legend(override.aes = list(size = .75)),
fill = guide_legend(override.aes = list(size = .75)),
size = "none", linetype = guide_legend(override.aes = list(size = .75))) +
coord_cartesian(ylim = c(0, 1), xlim = c(0,5))+
scale_y_continuous(breaks = seq(0,1,0.25)) +
scale_x_continuous(breaks = c(0,1,2,3,4,5))
2.2.4 Plot correlation and interval by item
nd <- crossing(domain= NA,
sei = 0.1,
item= unique(df$item),
event= NA,
formulation= NA,
age_dec_c= 0,
age_dec_c2= 0,
female_c= 0,
interval_val=seq(0,5, by = .1))
epred_draws_df <- nd %>%
add_epred_draws(fit_masc_m3, re_formula = NA)
# epred_draws_dom <- epred_draws_df %>%
# group_by(interval_val, health_subdomain) %>%
# mean_hdci(.epred,.width = c(.95,.8,.5)) %>%
# pivot_wider(names_from = .width, values_from = c(.lower,.upper))
epred_draws_df <- epred_draws_df %>%
mutate(item = str_to_title(as.character(item)))
ggplot(epred_draws_df) +
stat_lineribbon(alpha = 1/4, point_interval = "mean_hdci", aes(x = interval_val, y = .epred)) +
geom_point(data= df, aes(x=interval_val, y= cor_val), size = 1.5, color = "grey20", alpha = 0.7) +
# geom_line(data = epred_draws_agg,
# aes(x = time_diff_dec*10, y = .epred),
# color = "grey95",
# size = .5) +
# geom_line(data = epred_draws_dom,
# aes(x = time_diff_dec*10, y = .epred, linetype = domain_name),
# color = "#e07f00",
# linewidth = .25) +
# geom_text_repel(data = lbl_dot_df,
# aes(x = time_diff_dec*10, y = .epred, label = label),
# family = "Source Sans 3", size = 2.5,
# min.segment.length = 0,
# segment.color = "grey50",
# segment.size = .25,
# box.padding = 0.5,
# nudge_x = .5,
# nudge_y = c(.1, -.05)
# ) +
facet_wrap(.~str_to_title(item)) +
theme_minimal() +
labs(y = "Retest Correlation", x = "Retest Interval (Years)", color = "", linetype = "", fill = "",
title = "Behaviour") +
theme(strip.placement = "outside",
legend.position = "none", # c(0,0) bottom left, c(1,1) top-right.
legend.margin = margin(-.5,0,0,0, unit="cm"),
legend.spacing.y = unit(0.15, 'cm'),
legend.key.width = unit(1, "cm"),
legend.key.size = unit(.3, "cm"),
# plot.tag.position = c(0,.8),
legend.text = element_text(size = 10, color = "grey20"),
text = element_text(size = 12, color = "grey40"),
axis.text.y = element_text( vjust=seq(0,1, length.out = 5)),
axis.text.x = element_text( hjust=c(0,1)),
title = element_text(size = 12, color = "grey20"),
plot.tag = element_text( size = 12, face = "bold", color = "grey20"),
panel.spacing = unit(.5, "lines"),
plot.title = element_blank(),
panel.grid = element_blank(),
plot.title.position = "plot",
plot.margin = margin(b = 5, r = 5, l = 5),
panel.background = element_rect(color = "grey75", fill = NA, size = .4) ,
strip.text = element_text(face = "bold"),
strip.text.y = element_text(size = 12),
strip.text.x = element_text(size = 12))+
scale_fill_manual(values = c("#502a7b","#502a7c" , "#502a7a" )) +
guides(color = guide_legend(override.aes = list(size = .75)),
fill = guide_legend(override.aes = list(size = .75)),
size = "none", linetype = guide_legend(override.aes = list(size = .75))) +
coord_cartesian(ylim = c(0, 1), xlim = c(0,5)) +
scale_y_continuous(breaks = seq(0,1,0.25)) +
scale_x_continuous(breaks = c(0,1,2,3,4,5))
2.2.5 Plot correlation and interval by formulation
nd <- crossing(domain= NA,
sei = 0.1,
item= NA,
event= NA,
formulation= unique(df$formulation),
age_dec_c= 0,
age_dec_c2= 0,
female_c= 0,
interval_val=seq(0,5, by = .1))
epred_draws_df <- nd %>%
add_epred_draws(fit_masc_m3, re_formula = NA)
# epred_draws_dom <- epred_draws_df %>%
# group_by(interval_val, health_subdomain) %>%
# mean_hdci(.epred,.width = c(.95,.8,.5)) %>%
# pivot_wider(names_from = .width, values_from = c(.lower,.upper))
# Set order and capitalise first letter
epred_draws_df$formulation <- str_to_title(epred_draws_df$formulation)
epred_draws_df$formulation<- factor(epred_draws_df$formulation, levels = c("Risk","Concern", "Worry" ,"Threat"))
df_plot <- df %>%
mutate(formulation = str_to_title(formulation))
df_plot$formulation<- factor(df_plot$formulation, levels = c("Risk","Concern", "Worry" ,"Threat"))
# Plot
ggplot(epred_draws_df) +
stat_lineribbon(alpha = 1/4, point_interval = "mean_hdci", aes(x = interval_val, y = .epred)) +
geom_point(data= df_plot, aes(x=interval_val, y= cor_val), size = 1.5, color = "grey20", alpha = 0.7) +
# geom_line(data = epred_draws_agg,
# aes(x = time_diff_dec*10, y = .epred),
# color = "grey95",
# size = .5) +
# geom_line(data = epred_draws_dom,
# aes(x = time_diff_dec*10, y = .epred, linetype = domain_name),
# color = "#e07f00",
# linewidth = .25) +
# geom_text_repel(data = lbl_dot_df,
# aes(x = time_diff_dec*10, y = .epred, label = label),
# family = "Source Sans 3", size = 2.5,
# min.segment.length = 0,
# segment.color = "grey50",
# segment.size = .25,
# box.padding = 0.5,
# nudge_x = .5,
# nudge_y = c(.1, -.05)
# ) +
facet_wrap(.~formulation) +
theme_minimal() +
labs(y = "Retest Correlation", x = "Retest Interval (Years)", color = "", linetype = "", fill = "",
title = "Behaviour") +
theme(strip.placement = "outside",
legend.position = "none", # c(0,0) bottom left, c(1,1) top-right.
legend.margin = margin(-.5,0,0,0, unit="cm"),
legend.spacing.y = unit(0.15, 'cm'),
legend.key.width = unit(1, "cm"),
legend.key.size = unit(.3, "cm"),
# plot.tag.position = c(0,.8),
legend.text = element_text(size = 10, color = "grey20"),
text = element_text(size = 12, color = "grey40"),
axis.text.y = element_text( vjust=seq(0,1, length.out = 5)),
axis.text.x = element_text( hjust=c(0,1)),
title = element_text(size = 12, color = "grey20"),
plot.tag = element_text( size = 12, face = "bold", color = "grey20"),
panel.spacing = unit(.5, "lines"),
plot.title = element_blank(),
panel.grid = element_blank(),
plot.title.position = "plot",
plot.margin = margin(b = 5, r = 5, l = 5),
panel.background = element_rect(color = "grey75", fill = NA, size = .4),
strip.text = element_text(face = "bold"),
strip.text.y = element_text(size = 12),
strip.text.x = element_text(size = 12)) +
scale_fill_manual(values = c("#502a7b","#502a7c","#502a7a" ))+
guides(color = guide_legend(override.aes = list(size = .75)),
fill = guide_legend(override.aes = list(size = .75)),
size = "none", linetype = guide_legend(override.aes = list(size = .75))) +
coord_cartesian(ylim = c(0, 1), xlim = c(0,5))+
scale_y_continuous(breaks = seq(0,1,0.25)) +
scale_x_continuous(breaks = c(0,1,2,3,4,5))
2.3 Prediction Parameters
2.3.1 By Domain
pred_df_domain <- NULL
for (curr_nlpar in c("stabch","rel","change")) {
nd <- crossing(domain= unique(df$domain),
sei = 0.1,
item=NA,
event=NA,
formulation=NA,
age_dec_c=0,
age_dec_c2=0,
female_c= 0,
interval_val=0)
fit_nlpar_domain <- nd %>%
add_epred_draws(fit_masc_m3, nlpar = curr_nlpar, re_formula = NA)
fit_nlpar_domain <- fit_nlpar_domain %>%
group_by(domain) %>%
mean_hdci(.epred,.width = c(.95,.8,.5)) %>%
pivot_wider(names_from = .width, values_from = c(.lower,.upper)) %>%
mutate(nlpar = curr_nlpar,
estimate = .epred,
categ = "domain",
x = domain)%>%
select(categ, x, nlpar, estimate, dplyr::contains("er_"))
pred_df <- fit_nlpar_domain
pred_df_domain <- bind_rows(pred_df, pred_df_domain)
}
pred_df_domain <- pred_df_domain %>%
mutate(nlpar = case_when(nlpar == "rel" ~ "Reliability",
nlpar == "change" ~ "Change",
nlpar == "stabch" ~"Stab. Change"))2.3.2 By Item
pred_df_item <- NULL
for (curr_nlpar in c("stabch","rel","change")) {
nd <- crossing(domain= NA,
sei = 0.1,
item= unique(df$item),
event=NA,
formulation=NA,
age_dec_c=0,
age_dec_c2=0,
female_c= 0,
interval_val=0)
fit_nlpar_item <- nd %>%
add_epred_draws(fit_masc_m3, nlpar = curr_nlpar, re_formula = NA)
fit_nlpar_item <- fit_nlpar_item %>%
group_by(item) %>%
mean_hdci(.epred,.width = c(.95,.8,.5)) %>%
pivot_wider(names_from = .width, values_from = c(.lower,.upper)) %>%
mutate(nlpar = curr_nlpar,
estimate = .epred,
categ = "item",
x = item)%>%
select(categ, x, nlpar, estimate, dplyr::contains("er_"))
pred_df <- fit_nlpar_item
pred_df_item <- bind_rows(pred_df, pred_df_item)
}
pred_df_item <- pred_df_item %>%
mutate(nlpar = case_when(nlpar == "rel" ~ "Reliability",
nlpar == "change" ~ "Change",
nlpar == "stabch" ~"Stab. Change"))2.3.3 By Event
pred_df_event <- NULL
for (curr_nlpar in c("stabch","rel","change")) {
nd <- crossing(domain= NA,
sei = 0.1,
item= NA,
event=unique(df$event),
formulation=NA,
age_dec_c=0,
age_dec_c2=0,
female_c= 0,
interval_val=0)
fit_nlpar_event <- nd %>%
add_epred_draws(fit_masc_m3, nlpar = curr_nlpar, re_formula = NA)
fit_nlpar_event <- fit_nlpar_event %>%
group_by(event) %>%
mean_hdci(.epred,.width = c(.95,.8,.5)) %>%
pivot_wider(names_from = .width, values_from = c(.lower,.upper)) %>%
mutate(nlpar = curr_nlpar,
estimate = .epred,
categ = "event",
x = event)%>%
select(categ, x, nlpar, estimate, dplyr::contains("er_"))
pred_df <- fit_nlpar_event
pred_df_event <- bind_rows(pred_df, pred_df_event)
}
pred_df_event <- pred_df_event %>%
mutate(nlpar = case_when(nlpar == "rel" ~ "Reliability",
nlpar == "change" ~ "Change",
nlpar == "stabch" ~"Stab. Change"))2.3.4 By Formulation
pred_df_formulation <- NULL
for (curr_nlpar in c("stabch","rel","change")) {
nd <- crossing(domain= NA,
sei = 0.1,
item= NA,
event=NA,
formulation= unique(df$formulation),
age_dec_c=0,
age_dec_c2=0,
female_c= 0,
interval_val=0)
fit_nlpar_formulation <- nd %>%
add_epred_draws(fit_masc_m3, nlpar = curr_nlpar, re_formula = NA)
fit_nlpar_formulation <- fit_nlpar_formulation %>%
group_by(formulation) %>%
mean_hdci(.epred,.width = c(.95,.8,.5)) %>%
pivot_wider(names_from = .width, values_from = c(.lower,.upper)) %>%
mutate(nlpar = curr_nlpar,
estimate = .epred,
categ = "formulation",
x = formulation)%>%
select(categ, x, nlpar, estimate, dplyr::contains("er_"))
pred_df <- fit_nlpar_formulation
pred_df_formulation <- bind_rows(pred_df, pred_df_formulation)
}
pred_df_formulation <- pred_df_formulation %>%
mutate(nlpar = case_when(nlpar == "rel" ~ "Reliability",
nlpar == "change" ~ "Change",
nlpar == "stabch" ~"Stab. Change"))2.3.5 By Age
pred_df_age <- NULL
for (curr_nlpar in c("stabch","rel","change")) {
nd <- crossing(domain= NA,
sei = 0.1,
item= NA,
event=NA,
formulation=NA,
age_dec_c= c(-2,0,2,4),
age_dec_c2= c(4, 0, 4, 16),
female_c= 0,
interval_val=0)
fit_nlpar_age <- nd %>%
add_epred_draws(fit_masc_m3, nlpar = curr_nlpar, re_formula = NA)
fit_nlpar_age <- fit_nlpar_age %>%
group_by(age_dec_c) %>%
mean_hdci(.epred,.width = c(.95,.8,.5)) %>%
pivot_wider(names_from = .width, values_from = c(.lower,.upper)) %>%
mutate(nlpar = curr_nlpar,
estimate = .epred,
categ = "age",
x = case_when(age_dec_c== 4 ~ "80+", age_dec_c== 2 ~ "60-79", age_dec_c== 0 ~ "40-59", age_dec_c== -2 ~ "20-39"))%>%
select(categ, x, nlpar, estimate, dplyr::contains("er_"))
pred_df <- fit_nlpar_age
pred_df_age <- bind_rows(pred_df, pred_df_age)
}
pred_df_age <- pred_df_age %>%
mutate(nlpar = case_when(nlpar == "rel" ~ "Reliability",
nlpar == "change" ~ "Change",
nlpar == "stabch" ~"Stab. Change"))2.3.6 By Gender
pred_df_female <- NULL
for (curr_nlpar in c("stabch","rel","change")) {
nd <- crossing(domain= NA,
sei = 0.1,
item= NA,
event=NA,
formulation=NA,
age_dec_c= 0,
age_dec_c2= 0,
female_c= c(-0.5, 0.5),
interval_val=0)
fit_nlpar_female <- nd %>%
add_epred_draws(fit_masc_m3, nlpar = curr_nlpar, re_formula = NA)
fit_nlpar_female <- fit_nlpar_female %>%
group_by(female_c) %>%
mean_hdci(.epred,.width = c(.95,.8,.5)) %>%
pivot_wider(names_from = .width, values_from = c(.lower,.upper)) %>%
mutate(nlpar = curr_nlpar,
estimate = .epred,
categ = "gender",
x = ifelse(female_c== -0.5, "male", "female"))%>%
select(categ, x, nlpar, estimate, dplyr::contains("er_"))
pred_df <- fit_nlpar_female
pred_df_female <- bind_rows(pred_df, pred_df_female)
}
pred_df_female <- pred_df_female %>%
mutate(nlpar = case_when(nlpar == "rel" ~ "Reliability",
nlpar == "change" ~ "Change",
nlpar == "stabch" ~"Stab. Change"))2.3.7 Overall
pred_df_overall <- NULL
for (curr_nlpar in c("stabch","rel","change")) {
nd <- crossing(domain= NA,
sei = 0.1,
item= NA,
event=NA,
formulation=NA,
age_dec_c= 0,
age_dec_c2= 0,
female_c= 0,
interval_val=0)
fit_nlpar_overall <- nd %>%
add_epred_draws(fit_masc_m3, nlpar = curr_nlpar, re_formula = NA)
fit_nlpar_overall <- fit_nlpar_overall %>%
mean_hdci(.epred,.width = c(.95,.8,.5)) %>%
pivot_wider(names_from = .width, values_from = c(.lower,.upper)) %>%
mutate(nlpar = curr_nlpar,
estimate = .epred,
categ = "all",
x = "overall")%>%
select(categ, x, nlpar, estimate, dplyr::contains("er_"))
pred_df <- fit_nlpar_overall
pred_df_overall <- bind_rows(pred_df, pred_df_overall)
}
pred_df_overall <- pred_df_overall %>%
mutate(nlpar = case_when(nlpar == "rel" ~ "Reliability",
nlpar == "change" ~ "Change",
nlpar == "stabch" ~"Stab. Change"))
pred_df <- bind_rows(pred_df_domain, pred_df_item, pred_df_event, pred_df_formulation ,pred_df_age, pred_df_female, pred_df_overall)2.4 Parameters Table
library(gt)
# Selecting only the specified columns from pred_df
pred_df_selected <- pred_df[, c("categ", "x", "nlpar", "estimate", ".lower_0.95", ".upper_0.95")]
# Create a gt table
table_gt <- gt(data = pred_df_selected,
rowname_col = "x",
caption = "Prediction Results",
rownames_to_stub = TRUE)
# Formatting
table_gt <- table_gt %>%
tab_header(title = "Estimates and Confidence Intervals") %>%
fmt_number(columns = c("estimate", ".lower_0.95", ".upper_0.95"), decimals = 3) %>%
cols_align(align = "center")
table_gt| Estimates and Confidence Intervals | ||||||
| categ | x | nlpar | estimate | .lower_0.95 | .upper_0.95 | |
|---|---|---|---|---|---|---|
| 1 | domain | health | Change | 0.459 | 0.069 | 0.881 |
| 2 | domain | other | Change | 0.428 | 0.039 | 0.850 |
| 3 | domain | health | Reliability | 0.775 | 0.657 | 0.897 |
| 4 | domain | other | Reliability | 0.772 | 0.648 | 0.900 |
| 5 | domain | health | Stab. Change | 0.480 | 0.090 | 0.898 |
| 6 | domain | other | Stab. Change | 0.539 | 0.111 | 0.945 |
| 7 | item | multiple | Change | 0.425 | 0.054 | 0.872 |
| 8 | item | single | Change | 0.459 | 0.038 | 0.914 |
| 9 | item | multiple | Reliability | 0.745 | 0.624 | 0.864 |
| 10 | item | single | Reliability | 0.795 | 0.640 | 0.942 |
| 11 | item | multiple | Stab. Change | 0.408 | 0.033 | 0.862 |
| 12 | item | single | Stab. Change | 0.609 | 0.161 | 0.973 |
| 13 | event | none | Change | 0.394 | 0.042 | 0.826 |
| 14 | event | event | Change | 0.490 | 0.094 | 0.923 |
| 15 | event | none | Reliability | 0.838 | 0.739 | 0.932 |
| 16 | event | event | Reliability | 0.694 | 0.553 | 0.848 |
| 17 | event | none | Stab. Change | 0.461 | 0.055 | 0.886 |
| 18 | event | event | Stab. Change | 0.559 | 0.153 | 0.950 |
| 19 | formulation | risk | Change | 0.269 | 0.001 | 0.787 |
| 20 | formulation | concern | Change | 0.492 | 0.078 | 0.930 |
| 21 | formulation | threat | Change | 0.555 | 0.127 | 0.943 |
| 22 | formulation | worry | Change | 0.504 | 0.103 | 0.907 |
| 23 | formulation | risk | Reliability | 0.656 | 0.541 | 0.784 |
| 24 | formulation | concern | Reliability | 0.773 | 0.542 | 0.989 |
| 25 | formulation | threat | Reliability | 0.711 | 0.545 | 0.875 |
| 26 | formulation | worry | Reliability | 0.879 | 0.784 | 0.970 |
| 27 | formulation | risk | Stab. Change | 0.730 | 0.203 | 1.000 |
| 28 | formulation | concern | Stab. Change | 0.467 | 0.053 | 0.893 |
| 29 | formulation | threat | Stab. Change | 0.396 | 0.024 | 0.816 |
| 30 | formulation | worry | Stab. Change | 0.407 | 0.038 | 0.830 |
| 31 | age | 20-39 | Change | 0.583 | 0.011 | 1.000 |
| 32 | age | 40-59 | Change | 0.651 | 0.034 | 1.000 |
| 33 | age | 60-79 | Change | 0.710 | 0.029 | 1.000 |
| 34 | age | 80+ | Change | 0.737 | 0.011 | 1.000 |
| 35 | age | 20-39 | Reliability | 0.744 | 0.523 | 0.983 |
| 36 | age | 40-59 | Reliability | 0.766 | 0.574 | 1.000 |
| 37 | age | 60-79 | Reliability | 0.781 | 0.576 | 1.000 |
| 38 | age | 80+ | Reliability | 0.790 | 0.556 | 1.000 |
| 39 | age | 20-39 | Stab. Change | 0.496 | 0.000 | 0.998 |
| 40 | age | 40-59 | Stab. Change | 0.558 | 0.000 | 0.998 |
| 41 | age | 60-79 | Stab. Change | 0.614 | 0.000 | 0.999 |
| 42 | age | 80+ | Stab. Change | 0.646 | 0.000 | 1.000 |
| 43 | gender | male | Change | 0.435 | 0.083 | 0.830 |
| 44 | gender | female | Change | 0.448 | 0.093 | 0.835 |
| 45 | gender | male | Reliability | 0.752 | 0.603 | 0.893 |
| 46 | gender | female | Reliability | 0.793 | 0.672 | 0.911 |
| 47 | gender | male | Stab. Change | 0.523 | 0.141 | 0.892 |
| 48 | gender | female | Stab. Change | 0.496 | 0.103 | 0.857 |
| 49 | all | overall | Change | 0.439 | 0.114 | 0.805 |
| 50 | all | overall | Reliability | 0.775 | 0.667 | 0.894 |
| 51 | all | overall | Stab. Change | 0.510 | 0.159 | 0.861 |
2.5 Parameters Plot
pred_df <- pred_df %>%
mutate(categ = str_to_title(as.character(categ)),
x = str_to_title(as.character(x)))
pred_df$nlpar <- factor(pred_df$nlpar, levels = c("Reliability", "Change", "Stab. Change"))
pred_df$x <- factor(pred_df$x, levels = c("Overall", "None","Event", "Other", "Health", "Single", "Multiple", "Threat", "Concern", "Worry", "Risk", "80+","60-79", "40-59", "20-39", "Male", "Female"))
pred_df$categ <- factor(pred_df$categ, levels = c("All", "Event", "Domain", "Item", "Formulation" ,"Age", "Gender"))
p_nlpar <- pred_df %>% ggplot() +
geom_crossbar(aes(xmin = .lower_0.95, x = estimate,
xmax = .upper_0.95, y = x),
fill = "white", color = "NA",
linewidth = .15,width = 0.25, alpha = 1) +
geom_crossbar(aes(xmin = .lower_0.95, x = estimate,
xmax = .upper_0.95, y = x),
fill = "#502a7a", color = "NA",
linewidth = .15,width = 0.25, alpha = .3) +
geom_crossbar(aes(xmin = .lower_0.8, x = estimate,
xmax = .upper_0.8, y = x),
fill = "white", color = "NA",
linewidth = .15,width = 0.25, alpha = 1) +
geom_crossbar(aes(xmin = .lower_0.8, x = estimate,
xmax = .upper_0.8, y = x),
fill = "#502a7a",color = "NA",
linewidth = .15,width = 0.25, alpha = .6) +
geom_crossbar(aes(xmin = .lower_0.5, x = estimate,
xmax = .upper_0.5, y = x),
fill = "white",color = "NA",
linewidth = .15,width = 0.25, alpha = 1) +
geom_crossbar(aes(xmin = .lower_0.5, x = estimate,
xmax = .upper_0.5, y = x),
fill = "#502a7a",color = "NA",
linewidth = .15,width = 0.25, alpha = .9) +
geom_point(aes(x = estimate, y =x),
fill = "white", color = "grey20",
shape = 21,
stroke = .25,
size = 2) +
facet_grid(categ ~ nlpar, scales = "free_y", space = "free", switch = "y") +
# scale_x_continuous(limits = c(0,1), expand = c(0,0)) +
scale_x_continuous(limits = c(0, 1), expand = c(0, 0), breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) +
scale_y_discrete(position = "left") +
theme_minimal() +
geom_rect(data = subset(pred_df, x %in% c("Overall")),
fill = NA, colour = "black", size = .75, xmin = -Inf,xmax = Inf,
ymin = -Inf,ymax = Inf) +
theme(panel.grid = element_blank(),
legend.position = "none",
plot.title.position = "plot",
strip.placement = "outside",
#strip.text.y = element_blank(),
axis.title.x = element_text(size = 10, color = "grey20"),
plot.margin = margin(b = 10, t = 10, r = 15, l = 0),
panel.spacing.x = unit(.3, "cm"),
# plot.tag.position = c(0.025,.9),
# panel.grid.major.x = element_line(linewidth = .2, color = "grey80", linetype = "solid"),
panel.background = element_rect(linewidth = .25, color = "grey50", fill = "NA"),
#plot.background = element_rect(linewidth = .25, color = "grey40", fill = "NA"),
strip.text = element_text( face = "bold"),
strip.text.y = element_text(size = 10),
strip.text.x = element_text(size = 10),
plot.tag = element_text( size = 10, face = "bold", color = "grey20"),
plot.title = element_blank(),
title = element_text( face = "bold", size = 10),
axis.text.x = element_markdown( color = "black", size = 10, hjust=c(0,.5, 1)),
axis.text.y.left = element_markdown( angle = 0, hjust = 1, color = "grey20", size = 10), # hjust = c(0,.5,.5,.5,1)
text = element_text()) +
labs(y = "", x = "Parameter Estimate", title = "Propensity")
p_nlpar
ggsave(plot = p_nlpar,
filename = here::here("analysis", "masc_pred.jpeg"),
dpi = 300, width = 28, height = 30, units = "cm")